/*
CSE30332 Final Project - Iris Matcher
Group Oatmeal

Filter Program:
===================
This program runs a log-gabor filter over an iris image passed
as a command line argument.  After running the filter, it returns
one of two things (or both, depending on the options at command-line
parameters specified.  If the -o flag is specified, it will output
a bitvector; if the -v flag is specified, it will output a feature
vector.  The bitvector is computed by looking (pixel by pixel) at which
quadrant the result of the the IDFT (on the filtered image) is located in;
it either outputs 00, 01, 10, or 11 (or XX, if the bit is determine to be
occluded).  The feature vector simply takes the magnitude of the real and
imaginary parts of the filtered image (after running the DFT and IDFT) and
outputs this.

The log-gabor filter we ran is rather simple; we just run a 1-D DFT on the 
original image, run the log-gabor filter on this image, and then compute
the 1-D IDFT on it to give us an image with (hopefully) less noise. This is
essentially the same as the method Masek implemented in his study; the only
real difference are the values for f0 and sigma.  Masek used values of 36 and 18;
however, since our images were twice as large, we determined it would be most
effective to use double Masek's values. 
*/


#include <unistd.h>
#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>


void usage()
{
	printf("USAGE: ./filter -f filename [-o bitvector filename] [-v featurevector filename] [-d]\n");
	printf("\t-f the filename of the image to load (REQUIRED)\n");
	printf("\t-o the filename to write the bitvector out to\n");
	printf("\t-v the filename to write the featurevector out to\n");
	printf("\t-d prints debugging images to the screen\n");
}


int main (int argc, char * argv[])
{
        char *fname = 0x0;
        float f0 = 72.0;
        char *outname = 0x0;        
        char *fvname = 0x0;
        float sigma = 36.0;
        int debug = 0;
        char c;
        FILE *out2b, *outfv;
        int i,j,a,b;
        int mean;
        char bf = 0;
	char vf = 0;
	double dev;
	
        //handle the various arguments
        while ((c = getopt( argc, argv, ":f:o:v:d" )) != -1)
        {
                switch(c)
                {
                        case 'f':
                                fname = optarg;
                                break;
                        case 'o':
                                outname = optarg;
				bf++;
                                break;
                        case 'v':
                                fvname = optarg;
				vf++;
				break;
                        case 'd':
                                debug++;
                                break;
                        case '?':
                        default:
				usage();
                                fprintf(stderr, "Error; Need filename, method\n");
                }
        }

        if (fname == 0x0)
        {
		usage();
                fprintf(stderr, "Invalid filename\n");
                exit(-1);
        }



	if (bf)
        	out2b = fopen(outname, "w");

        if (vf)
                outfv = fopen(fvname, "w");
        
        IplImage *src = cvLoadImage(fname,0);
        
	if (!src)
	{
		printf("Error loading image\n");
		exit(-1);
	}
	
	if (debug) 
        {
                cvNamedWindow("Source", 1);
                cvShowImage("Source", src);
                cvWaitKey(0);
        }

	/*allocate the memory for the DFT images*/
        IplImage *inReal = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_32F, 1);
        cvConvert(src, inReal);
        IplImage *dftIn = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_32F, 2);
        cvMerge(inReal, NULL, NULL, NULL, dftIn);
        IplImage *dftOut = cvCreateImage(cvGetSize(dftIn), IPL_DEPTH_32F, 2);
	/*perform the actual DFT*/
        cvDFT(dftIn, dftOut, CV_DXT_FORWARD|CV_DXT_ROWS, 0);

	/*show the results of the DFT to screen*/
        if (debug)
        {
                IplImage *outRe = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
                IplImage *outIm = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
                cvSplit(dftOut, outRe, outIm, NULL, NULL);
                cvNamedWindow("After DFT Real", 1);
                cvShowImage("After DFT Real", outRe);
                cvNamedWindow("After DFT Im", 1);
                cvShowImage("After DFT Im", outIm);
                cvWaitKey(0);
        }

	/*Set up the log-gabor filter*/
        float d = 2.0*log(sigma/f0) * log(sigma/f0);
        int w = src->width;
        int h = src->height;
        IplImage *fOut = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 2);
        
	/*Perform the filtering operation on each individual pixel in the image*/
        for (j = 0; j < w; j++)
        {
                float f = abs(1.0*j-w/2.0);
                float mod = exp(-1.0*log(f/f0)*log(f/f0)/d);
                if (isnan(mod)) mod = 0.0;
                for (i = 0; i < h; i++)
                {
                        CvScalar s = cvGet2D(dftOut,i,j);
                        s.val[0] *= mod;
                        s.val[1] *= mod;
                        cvSet2D(fOut,i,j,s);
                }
        }

	/*show the results of the filter to the screen*/
        if (debug)
        {
                IplImage *foutRe = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
                IplImage *foutIm = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
                cvSplit(fOut, foutRe, foutIm, NULL, NULL);
                cvNamedWindow("After Filter Real", 1);
                cvShowImage("After Filter Real", foutRe);
                cvNamedWindow("After Filter Im", 1);
                cvShowImage("After Filter Im", foutIm);
                cvWaitKey(0);
        }


        /*allocate the necessary memory and perform the inverse DFT*/ 
        IplImage *idftOut = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 2);
        cvDFT(fOut, idftOut, CV_DXT_INV_SCALE|CV_DXT_ROWS,0);
        IplImage *ioutRe = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
        IplImage *ioutIm = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_32F, 1);
        cvSplit(idftOut, ioutRe, ioutIm, NULL, NULL);

	/*output the results of the IDFT to screen*/
        if (debug)
        {
                IplImage *re = cvCreateImage(cvGetSize(dftOut), IPL_DEPTH_8U, 1);
                cvConvert(ioutRe, re);
                cvNamedWindow("After IDFT Real", 1);
                cvShowImage("After IDFT Real", re);
                cvNamedWindow("After IDFT Im", 1);
                cvShowImage("After IDFT Im", ioutIm);
                cvWaitKey(0);
        }

	/*allocate the space for the final images*/
        IplImage *acttemp = cvCreateImage(cvGetSize(ioutRe), IPL_DEPTH_32F, 1);
        IplImage *act = cvCreateImage(cvGetSize(ioutIm), IPL_DEPTH_8U, 1);

        double val;
	double origval;
        double yi,yr;

	/*write the iris codes or mask values*/
        for (j = 0; j < w; j++)
        {
                for (i = 0; i < h-4; i++)
                { 
			/*compute the magnitude of the result of the IDFT (since it is a complex number)*/
                        yi = cvGetReal2D(ioutIm,i,j);
                        yr = cvGetReal2D(ioutRe,i,j);
			val = yi*yi+yr*yr;
                        val = sqrt(val);
                        
			/*if the value in the original is above or below a threshold, output masking bit*/
			origval = cvGetReal2D(src,i,j);
			if (origval >= 240)
			{
				if (bf) fprintf(out2b, "XX");
                                if (vf) fprintf(outfv, "0\n");
                        }
                        else if (origval <= 6)
			{
                                if (bf) fprintf(out2b, "XX");	
                                if (vf) fprintf(outfv, "0\n");
                        }
                        else
			{
				if (bf)
				{
	                        	/*write 2bit iris codes*/
        	        		if (yi > 0)
                	        	        fprintf(out2b, "0");
                        		else    
	                        	        fprintf(out2b, "1");
        	                	if (yr > 0)
                	        	        fprintf(out2b, "0");
                        		else
                                		fprintf(out2b, "1");
				}
				if (vf) fprintf(outfv, "%f\n",val);
			}
			/*set the image for debugging purposes*/
                        cvSetReal2D(acttemp, i,j,val);
                        val = 0;
                }
        }       
        
        if (vf) fclose(outfv);
        if (bf) fclose(out2b);
        
	/*for debugging purposes, normalize the image and output it to screen*/
        if (debug)
        {
        	cvNormalize(acttemp, acttemp, 255.0, 0, CV_MINMAX, NULL);
        	cvConvert(acttemp, act);
                cvNamedWindow("Final", 1);
                cvShowImage("Final", act);
                cvWaitKey(0);
        }
        return 0;
}
